\documentclass[a4paper]{article}
\usepackage[affil-it]{authblk}
\usepackage{indentfirst}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{geometry}
\usepackage{xeCJK}
\usepackage{enumitem}
\geometry{margin=1.5cm, vmargin={0pt,1cm}}
\setlength{\topmargin}{-1cm}
\setlength{\paperheight}{29.7cm}
\setlength{\textheight}{25.3cm}

\begin{document}
% ==================================================
\title{Numerical Analysis homework \# 6}

\author{周川迪 Zhou Chuandi 3220101409}
\affil{强基数学2201}

\date{\today}

\maketitle

\begin{abstract}
6.6.1 Theoretical questions I-IV from \textit{handouts NUMPDEs-2024-08-17}.
\end{abstract}

\section*{I. Simpson's rule}

\subsection*{(a)}
Simpson's rule on \(\int_{-1}^1 y(t) \, \mathrm{d}t\) with \(\rho=1\) is \(I^S(y) = \frac{1}{3} (y(-1) + 4y(0) + y(1)) \).
\begin{align*}
  &\ \int_{-1}^1 p_3(y;-1,0,0,-1; t) \, \mathrm{d}t \\
  =&\ \int_{-1}^1 \left[ y(-1) + y[-1,0](t+1) + y[-1,0,0] t(t+1) + y[-1,0,0,1](t+1)^2 t \right] \, \mathrm{d}t \\
  =&\ 2y(-1) + 2(y(0) - y(-1)) + \frac{2}{3} \left( y'(0) - y(0) + y(-1) \right) + \frac{2}{3} \cdot \frac{y'(0) - y(0) + y(-1) - y(1) + y(0) + y'(0)}{2}\\
  =&\ \frac{1}{3} (y(-1) + 4y(0) + y(1)) 
\end{align*}

Thus, the Simpson's rule on \([-1,1]\) with \(\rho=1\) can be obtained by \(\int_{-1}^1 y(t) \, \mathrm{d}t = \int_{-1}^1 p_3(y;-1,0,0,1;t) \, \mathrm{d}t + E^S(y)\).


\subsection*{(b)}
\[
E^S(t) = \int_{-1}^1 y(t) \, \mathrm{d}t - I^S(y) = \int_{-1}^1 y(t) \, \mathrm{d}t - \frac{1}{3} (y(-1) + 4y(0) + y(1)) = \int_{-1}^1 (y(t)-p_3(y;-1,0,0,1;t)) \, \mathrm{d}t
\]

From \textbf{Theorem 2.37}, \[y(t)-p_3(y;-1,0,0,1;t)=\frac{1}{4!}y^{(4)}(\xi)(t+1)t^2(t-1), \xi \in (-1,1)\] 

Thus, we have (for some \(\xi \in (-1,1)\))
\[E^S(y) = \frac{1}{24} y^{(4)}(\xi) \int_{-1}^1 (t+1)t^2(t-1) \, \mathrm{d}t = -\frac{1}{24} y^{(4)}(\xi) \cdot \frac{4}{15} = -\frac{1}{90} y^{(4)}(\xi)\]


\subsection*{(c)}
\subsubsection*{Proof of composite Simpson's rule: \(I_n^S(y)=\frac{h}{3} \left( y(x_0) + 4y(x_1) + 2y(x_2) + 4y(x_3) + \cdots + 4y(x_{n-2}) + y(x_n) \right)\)}

For \(h=\frac{b-a}{n}, x_k=a+kh\), apply the result in (a) to the subintervals \([x_{2k}, x_{2k+2}]\) and for \(k=0,1,\ldots,n-2\), and sum them up:
\begin{align*}
  I_n^S(y) &= \int_{x_0}^{x_2} p_3\left(y; x_0, x_1, x_1, x_2; t\right) \, \mathrm{d}t + \int_{x_2}^{x_4} p_3\left(y; x_2, x_3, x_3, x_4; t\right) \, \mathrm{d}t + \cdots + \int_{x_{n-2}}^{x_n} p_3\left(y; x_{n-2}, x_{n-1}, x_{n-1}, x_n; t\right) \, \mathrm{d}t \\
  &=\frac{h}{3} \left( y(x_0) + 4y(x_1) + y(x_2) \right) + \frac{h}{3} \left( y(x_2) + 4y(x_3) + 4y(x_4) \right) + \cdots + \frac{h}{3} \left( y(x_{n-2}) + 4y(x_{n-1}) + y(x_n) \right) \\
  &= \frac{h}{3} \left( y(x_0) + 4y(x_1) + 2y(x_2) + 4y(x_3) + \cdots + 4y(x_{n-2}) + y(x_n) \right).
\end{align*}

\subsubsection*{Proof of the error estimation: \(E_n^S(y) = -\frac{b-a}{180} h^4 y^{(4)}(\xi)\)}

Apply the result in (b) to each subinterval and sum them up:

\begin{align*}
  E_n^S(y) 
  &= \sum_{k=0}^{\frac{n}{2}-1}\int_{x_{2k}}^{x_{2k+2}}(y(t)- p_3(y; x_{2k}, x_{2k+1}, x_{2k+1}, x_{2k+2}; t))\, \mathrm{d}t \\
  &= \sum_{k=0}^{\frac{n}{2}-1} \int_{x_{2k}}^{x_{2k+2}} \frac{1}{4!} y^{(4)}(\xi_k) (t-x_{2k})(t-x_{2k+1})^2(t-x_{2k+2}) \, \mathrm{d}t \\
  &= \sum_{k=0}^{\frac{n}{2}-1} \frac{1}{24} y^{(4)}(\xi_k) h \int_{-h}^{h} (s+1)s^2(s-1) \, \mathrm{d}s \\
  &= \sum_{k=0}^{\frac{n}{2}-1} -\frac{y^{(4)}(\xi_k)}{90} \cdot h^5 \\
  &= -\frac{y^{(4)}(\xi)}{90} \cdot h^5 \cdot \frac{n}{2} \\
  &= -\frac{b-a}{180} h^4 y^{(4)}(\xi) \\
\end{align*}
  



\section*{II. Estimate the number of subintervals required}

We have \(f(x)=e^{-x^{2}}\) on \([0,1]\), \(h=\frac{b-a}{n}=\frac{1}{n}\).

\subsection*{(a)}
by the composite trapezoidal rule, \[E_n^T(f)= -\frac{b-a}{12} h^2 f''(\xi) = -\frac{1}{12n^2} f''(\xi)\].

We have \(f''(x)=e^{-x^2} (-2 + 4 x^2)\), so \(|f''(x)| \leq 2\) on \([0,1]\). Therefore, 
\[|E_n^T(f)| \leq \frac{1}{6n^2}\]

We want \(|E_n^T(f)| \leq 0.5 \times 10^{-6}\), so we need \(n \geq 578\).


\subsection*{(b)}
by the composite Simpson's rule, \[E_n^S(f) = -\frac{b-a}{180} h^4 f^{(4)}(\xi) = -\frac{1}{180n^4} f^{(4)}(\xi)\].

We have \(f^{(4)}(x)=e^{-x^2} (12 - 48 x^2 + 16 x^4)\), so \(|f^{(4)}(x)| \leq 12\) on \([0,1]\). Therefore,
\[|E_n^S(f)| \leq \frac{1}{15n^4}\]

We want \(|E_n^S(f)| \leq 0.5 \times 10^{-6}\), so we need \(n \geq 20\).

(Note: \(-2 + 4 x^2\) and \(12 - 48 x^2 + 16 x^4\) are ``opening upwards" on \(x \in [0,1]\), so we only need to consider their values at 0 and 1 to determine the upper bound.)






\section*{III. Gauss-Laguerre quadrature formula}

(hint: \(\int_{0}^{+\infty} t^{m} e^{-t} \, \mathrm{d}t = m!\))

\subsection*{(a)}
To construct \(\pi_{2}(t)=t^2+a t+b\) such that for \(\rho(t)=e^{-t}\), \(\forall p\in \mathbb{P}_1,\quad\int_0^{+\infty} p(t)\pi_2(t)\rho(t)dt=0\), 

It is sufficient to have \[\left\{
\begin{array}{ll}
  \int_{0}^{+\infty} t e^{-t} \pi_{2}(t) \, \mathrm{d}t = 0 \\
  \int_{0}^{+\infty} e^{-t} \pi_{2}(t) \, \mathrm{d}t = 0
\end{array}\right. \Rightarrow \left\{
\begin{array}{ll}
  \int_{0}^{+\infty} (t^3+a t^2+ bt) e^{-t}\, \mathrm{d}t = 0 \\
  \int_{0}^{+\infty} (t^2+a t+ b) e^{-t}\, \mathrm{d}t = 0 \\
\end{array}\right.\]

For polynomial \(P(t)\), \begin{align*}
  \int_{0}^{+\infty} P(t) e^{-t} \, \mathrm{d}t 
  &= -P(t) e^{-t} \big|_{0}^{+\infty} + \int_{0}^{+\infty} e^{-t} \, \mathrm{d}P(t) \\
  &= P(0) + \int_{0}^{+\infty} P'(t) e^{-t} \, \mathrm{d}t \\
  &= P(0) + P'(0) + P''(0) + \cdots
\end{align*}

Thus,\[\left\{
\begin{array}{ll}
  0+b+2a+6 = 0 \\
  b+a+2 = 0
\end{array}\right. \Rightarrow \left\{
\begin{array}{ll}
  a=-4 \\
  b=2
\end{array}\right.\Rightarrow \pi_2(t) = t^2-4t+2\]

\subsection*{(b)}
We want \(\pi_2(t)=(t-t_1)(t-t_2)\), so we have \(t_1=2+\sqrt{2}, t_2=2-\sqrt{2}\).

Then, to satisfy \(\forall f \in \mathbb{P}_1,\ \int_0^{+\infty} f(t) e^{-t} \, \mathrm{d}t = w_1 f(t_1) + w_2 f(t_2)\),

\[\left\{
\begin{array}{ll}
  f=1 \Rightarrow \int_0^{+\infty} e^{-t} \, \mathrm{d}t = w_1 + w_2 \\
  f=t \Rightarrow \int_0^{+\infty} t e^{-t} \, \mathrm{d}t = w_1 t_1 + w_2 t_2
\end{array}\right. \Rightarrow \left\{
\begin{array}{ll}
  w_1 + w_2 = 1 \\
  (2+\sqrt{2})w_1 + (2-\sqrt{2})w_2 = 1
\end{array}\right. \Rightarrow \left\{
\begin{array}{ll}
  w_1 = \frac{2+\sqrt{2}}{4} \\
  w_2 = \frac{2-\sqrt{2}}{4}
\end{array}\right.\]

Therefore, we have \[\int_0^{+\infty}f(t)e^{-t} \, \mathrm{d}t = \frac{2+\sqrt{2}}{4}f(2-\sqrt{2}) + \frac{2-\sqrt{2}}{4}f(2+\sqrt{2}) + E_2(f)\]

By \textbf{Theorem 6.33}, there exists some \(\tau \in (0, +\infty)\) such that \[E_2(f) = \frac{f^{(4)}(\tau)}{4!} \int_0^{+\infty} \rho(t) \pi_2^2(t) \, \mathrm{d}t = \frac{f^{(4)}(\tau)}{4!} \int_0^{+\infty} e^{-t} (t^2-4t+2)^2 \, \mathrm{d}t=\frac{f^{(4)}(\tau)}{6}\]



\subsection*{(c)}
(hint: use the exact value \(I = \int_0^{+\infty} \frac{1}{1+t} e^{-t} \, \mathrm{d}t= 0.596347361\cdots\))

From (b), we have that the approximation of $I$ is \[I_2 = \frac{2+\sqrt{2}}{4}\cdot \frac{1}{1+(2-\sqrt{2})} + \frac{2-\sqrt{2}}{4}\cdot \frac{1}{1+(2+\sqrt{2})} = \frac{4}{7} \approx 0.571428571\]

Since \(f(t)=\frac{1}{1+t}\), \(f^{(4)}(t)=\frac{24}{(1+t)^5}\), the remainder is \(E_2(f) = \frac{f^{(4)}(\tau)}{6} = \frac{4}{(1+\tau)^5}\).

With the exact value, the error \(E \approx 0.596347361 - 0.571428571 = 0.02491879\).

Let $E = E_2(f) = \frac{4}{(1+\tau)^5}$, we have \(\tau \approx 1.76126\)





\section*{IV. Remainder of Gauss formulas}



\subsection*{(a)}

For the elementary Lagrange polynomial \(\ell_m\), we have
\[\ell_m(x) = \prod_{j\neq m} \frac{x-x_j}{x_m-x_j},\ \ell_m(x_i) = \delta_{m,i}, \ell_m'(x_m) = \sum_{j \neq m} \frac{1}{x_m-x_j}\]

For \(p(t)=\sum_{m=1}^n [h_m(t) f_m + q_m(t) f_m'], \quad  h_m(t)=(a_m+b_mt)\ell_m^2(t), \quad q_m(t)=(c_m+d_mt)\ell_m^2(t)\),

we have \( h'(t) = b_m\ell_m^2(t) + (a_m+b_mt) 2\ell_m(t) \ell_m'(t), \quad q'(t) = d_m\ell_m^2(t) + (c_m+d_mt) 2\ell_m(t) \ell_m'(t). \)

Then,

\begin{align*}
  p(x_m) &= \sum_{k=1}^n \left[ h_k(x_m) f_k + q_k(x_m) f_k' \right] \\
  &= \sum_{k=1}^n \left[ (a_k + b_k x_m) f_k + (c_k + d_k x_m) f_k' \right] \ell_k^2(x_m) \\
  &= (a_m + b_m x_m) f_m + (c_m + d_m x_m) f_m'\\
  p'(x_m) &= \sum_{k=1}^n \left[ h_k'(x_m) f_k + q_k'(x_m) f_k' \right] \\
  &= \sum_{k=1}^n \left[ (b_k f_k + d_k f_k') \ell_k(x_m) + \left[ (a_k + b_k x_m) f_k + (c_k + d_k x_m) f_k' \right] 2 \ell_k'(x_m) \right] \ell_k(x_m) \\
  &= (b_m f_m + d_m f_m') + 2 \left[ (a_m + b_m x_m) f_m + (c_m + d_m x_m) f_m' \right] \ell_m'(x_m) \\
\end{align*}

We aim to find \(p \in P_{2n-1}\) such that \(p(x_m)=f_m,\ p'(x_m)=f_m'\). So,

\[
\left\{\begin{array}{ll}
  a_m + b_m x_m = 1 \\
  c_m + d_m x_m = 0 \\
  2(a_m + x_m b_m) \ell_m'(x_m) + b_m = 0 \\
  2(c_m + x_m d_m) \ell_m'(x_m) + d_m = 1
\end{array}\right. \Rightarrow \left\{\begin{array}{ll}
  a_m = 1 + 2x_m \ell_m'(x_m) = 1 + 2x_m \sum_{i \neq m} \frac{1}{x_m - x_i} \\
  b_m = -2 \ell_m'(x_m) = -2 \sum_{i \neq m} \frac{1}{x_m - x_i} \\
  c_m = -x_m \\
  d_m = 1
\end{array}\right. 
\]




\subsection*{(b)}
To obtain \( I_n(f) = \sum_{k=1}^n [w_kf(x_k)+\mu_kf'(x_k)]\), we just take \begin{align*}
  w_k &= \int_a^b h_k(t) \rho(t)\, \mathrm{d}t = \int_a^b (a_k+b_kt)\ell_k^2(t) \rho(t)\, \mathrm{d}t = \int_a^b \left[ 1 + 2(x_k - t) \ell_k'(x_k) \right] \ell_k^2(t) \rho(t)\, \mathrm{d}t \\
  \mu_k &= \int_a^b q_k(t) \rho(t)\, \mathrm{d}t = \int_a^b (c_k+d_kt)\ell_k^2(t) = \int_a^b (t - x_k) \ell_k^2(t) \rho(t)\, \mathrm{d}t
\end{align*}

Then, \begin{align*}
  \sum_{k=1}^n [w_kf(x_k)+\mu_kf'(x_k)] 
  &= \sum_{k=1}^n \int_a^b h_k(t)f(x_k) \rho(t)\, \mathrm{d}t + \sum_{k=1}^n \int_a^b q_k(t)f'(x_k) \rho(t)\, \mathrm{d}t \\
  &= \int_a^b \sum_{k=1}^n ( h_k(t)f(x_k) + q_k(t)f'(x_k)) \rho(t)\, \mathrm{d}t \\
  &= \int_a^b p(t) \rho(t)\, \mathrm{d}t
\end{align*}

Thus, \(\int_a^b f(t) \rho(t)\, \mathrm{d}t = I_n(f) + E_n(f) = \int_a^b p(t) \rho(t)\, \mathrm{d}t + E_n(f) \), where \( p \) is the Hermite interpolation of \( f \) of degree \( 2n-1 \) by (a).

For all \( p \in \mathbb{P}_{2n-1} \), any Hermite interpolation of \( p \) of degree \( 2n-1 \) is also \( p \), so \( E_n(p) = 0 \) for this rule.



\subsection*{(c)}
Since \(\mu_k = \int_a^b (t - x_k) l_k^2(t) \rho(t) \mathrm{d}t\), to make $\mu_k = 0$ for each $k = 1, 2, \dots, n$, we need 
\[0 = \int_a^b (t - x_k) l_k^2(t) \rho(t) \mathrm{d}t = \int_a^b \frac{\prod_{j=1}^{n}(t - x_j)}{\prod_{j \neq k}(x_k - x_j)} l_k(t) \rho(t) \mathrm{d}t =
\frac{1}{\prod_{j \neq k}(x_k - x_j)} \int_a^b v_n(t) l_k(t) \rho(t) \mathrm{d}t\]

i.e. \(\int_a^b v_n(t)l_k(t) \rho(t) \mathrm{d}t = 0\) for all \(k = 1, 2, \dots, n\)

Since $\{l_1(t), l_2(t), \dots, l_n(t)\}$ is a basis of $\mathbb{P}_{n-1}$, the above condition is equivalent to

$$
\forall p \in \mathbb{P}_{n-1}, \quad \int_a^b v_n(t)p(t) \rho(t) \mathrm{d}t = 0,
$$

which is the same as the result in \textbf{Theorem 6.24}.




\end{document}